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ABSTRACT 

We study the vertical settling of solid bodies in a turbulent protoplanetary disc. We 
consider the situation when the coupling to the gas is weak or equivalently when the 
particle stopping time r s < due to friction with the gas is long compared to the orbital 
^ • timescale f2 _1 . An analytical model, which takes into account the stochastic nature of 

| the sedimentation process using a Fokker-Planck equation for the particle distribution 

1"**- ■ function in phase space, is used to obtain the vertical scale height of the solid layer as 

' a function of the vertical component of the turbulent gas velocity correlation function 

and the particle stopping time. This is found to be of the same form as the relation 
obtained for strongly coupled particles in previous work. 

We compare the predictions of this model with results obtained from local shearing 
box MHD simulations of solid particles embedded in a vertically stratified disc in 
which there is turbulence driven by the MRI. We find that the ratio of the dust 
Q-f disc thickness to the gas disc thickness satifies Hd/H = 0.08(f2T s t) _1 / 2 , which is 

q in very good agreement with the analytical model. By discussing the conditions for 

' gravitational instability in the outer regions of protoplanetary discs in which there is 

^ \ a similar level of turbulence, we find that bodies in the size range 50 to 600 metres 

can aggregate to form Kuiper belt-like objects with characteristic radii ranging from 
tens to hundreds of kilometres. 
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1 INTRODUCTION 

Dust is the constituent of planetesimals, that are believed to 
lead to the formation of planets in our own and other solar 
systems. Observational evidence for growth of dust particles 
in protoplanetary discs comes from near and mid-infrared 
imaging, mid-infrared spectrometry and millimetre inter- 
ferometry llNatta et alJl200d iRodmann et af]l200dh They 
point to the presence of grains of millimetre and even cen- 
timeter size in discs around clasical T Tauri stars. 

In laminar discs, it is well known that gas drag causes 
the particles to sediment towards the disc midplane. In the 
absence of turbulence, there is nothing to oppose sedimen- 
tation and a very thin dust sub-layer can form. Thus there 
is a possibility that gravitational instability ocurs giving 
rise to gravitationally b ound clumps, eventually leading to 
planetesimal formation dGoldreich fc Wardlll973() . But col- 
lective effects arising through the interaction of the optically 
thick sedimenting dust layer with the gas may lead to the 
generation of turbu lence and inhibit gravitational instabil- 
ity (see for exan njtelWgid^mjChfl^ alJll993t 
lYoudin fc ShulhoOllOomez fc Ostrikerll200.^~ 

On the other hand, provided an adequate degree of ioni- 



sation can be maintained, accretion flows in which the dom- 
inant motion is Keplerian rotation, such as occur in pro- 
toplanetary discs, are not laminar. They have been shown 
to develop turbu lence as a result of the magnetorotational 
instability (MRI: lHawlev et al]ll995F) . The gaseous velocity 
fluctuations associated with the turbulence affect the spatial 
distribution of particles. Small particles are strongly coupled 
to the fluid and essentially follow the gas. 

T he ra dial diff usion of a p assive contaminant was stud- 
ied bv ICarballido et alJ J2005I) . who found a diffusion co- 
efficient approximately 10 times smaller than the effective 
disc kinematic viscosity associated with angular momen - 
tum transport. As shown recently bv lJohansen et all |2006), 
this ratio, found to be o f order unity in other simulations 
Jjohansen fc Klahjl200^ . is dependent on the topology of 
the initial magnetic field and the amount of conserved mag- 
netic flux present, which determines the level of the resulting 
turbulence. They note that the form of the particle diffu- 
sion induced by the turbulence is apparently determined by 
this level, or e quivalently, the am ount of angular momen- 
tum transport. iTurner et alJ i2006t) also studied the vertical 
spreading of a trace species in the upper layers of a stratified 
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disc. They found radial spreading to be faster than vertical 
spreading throughout the vertical extent of the disc. 

Larger particles, for which the coupling is weaker, are 
found to start to sediment despite the presence of turbu- 
lence. For a surface density £ = 200 g cm -2 being charac- 
teristic of the minimum mass solar nebula at 5.2 AU, sed- 
imentation was found to begin for particle sizes between a 
centimetre and a metre in a disc model with MRI driven 
turbulence iFromang fc Papaloizoull2006t hereafter FP). FP 
modelled this settling process by means of an advection- 
diffusion equation for the dust density. They expressed the 
diffusion coefficient in terms of the turbulent gas velocity 
correlation function. They found that turbulence stirs the 
dust sub-layer significantly, up to approximately 20% of the 
disc scale height for particles of about 10 cm in radius at 5.2 
AU. Thus dust particles of this size range could not be the 
constituents of a gravitationally unstable sub-disc. Accord- 
ingly, on a scale of 5 AU in a minimum mass solar nebula, the 
formation of planetesimals is hindered until a much larger 
particle size range is reached when the MRI operates. 

The above studies modelled the dust component as a 
second fluid. However, to study the effect of MHD turbu- 
lence on solids of larger size such that the stopping time due 
to gas drag significantly exceeds the inverse orbital angular 
frequency, it is necessary to model dust grains as discrete 
particles. Johansen et al. (2006) have performed local three- 
dimensional simulations of a non-stratified turbulent proto- 
planetary disc and find that individual particles tend to con- 
centrate by a factor of up to 100 due to radial concentration 
in vortices. They suggest that such dust agglomerations can 
potentially become gravitationally unstable. However, their 
study did not address the effect of the vertical sedimentation 
of solid bodies. 

It is the purpose of this paper to perform a study of the 
stochastic vertical settling of solids that is applicable when 
the stopping time is much larger than the orbital period. 
Firstly we develop an analytical model based on solving a 
one dimensional Fokkcr-Planck equation. The derivation of 
this equation and its solution are presented in Section 2. 
This gives an expression for the dust sub-disc thickness as a 
function of particle size. In Section 3, we go on to compare 
these results, and others found in previous work applicable 
to strongly coupled small particles, with those found from 
simulations of solid particles embedded in a turbulent pro- 
toplanetary disc. Simulations have been performed which 
treat the dust dynamics us ing the Epstein and Stokes drag 
laws IWeidenschillins|ll977l) . Note that these include locally 
induced radial motion, but not that induced by a global ra- 
dial pressure gardient, so that in the former context, the 
vertical settling that is found will occur regardless of radial 
concentration in vortices as long as the particles do not col- 
lide. Finally, in Section 4 we discuss the implications of our 
results for the formation of planetesimals and the possibility 
of gravitational instability in the outer regions of protoplan- 
etary discs. 



2 THE MODEL 

In a protoplanetary disc, solid bodies evolve under the com- 
bined influence of gas drag and MHD turbulence. By damp- 
ing motion relative to the gas, on average the former results 



in settling towards the equatorial plane of the disc. Simulta- 
neously, the latter excites random motions that oppose the 
settling process. The ensemble of solid bodies will eventu- 
ally reach a steady state with an associated vertical semi- 
thickness Hd- Below we present a model to determine this 
equilibrium state which is applicable when the coupling of 
the particle to the gas is weak and such that the stopping 
time is very much longer than the orbital period. 

2.1 Definitions and notation 

The interaction between the gas and the solid bodies (dust, 
boulders, planetesimals) occurs through a drag force Fn- 
This can be written as 



F D = m p - 



(1) 



where m p is the mass of a particle, v its velocity and u is 
the gas velocity. 

The characteristic time required to equilibrate gas and 
particle velocities, or stopping time is r st . This depends on 
the ratio of the size of the particle, expressed through the 
radius a they would have if they were sp herical, and the 
mean free path of the gas molecules, A iWeidenschillinel 
Il977t ICuzzi et al.lll993l) . When A > (4/9)o, we are in the 
Epstein regime for which 



JEp) 



p C s 



(2) 



where p s is the solid mass density and c s is the gas speed of 
sound. 

When A < (4/9)a, we are in the Stokes regime, and 



(St) 

r st is defined through 



_(st) _ 8 p s 



3 p Cd\u ■ 



(3) 



where Cd is a dimensionless coefficient that depends on the 
Reynolds number of the flow lZ e , tending to a constant value 
of order unity for lZ e — > oo. 

2.2 Diffusion in velocity— space 

Since we are interested in the problem of vertical settling, we 
restrict the analysis to consider only the vertical direction 
z, though we shall consider relevant results of simulations 
relating to the horizontal directions below. 

The vertical component of the equation of motion for a 
single particle moving under the gravitational acceleration 
due to the central star and gas drag can be written as 

dv > = -tf z+ ?>=l* (4 ) 



dt 



where u z is the vertical component of the gas velocity which 
is fluctuating stochastically because it is turbulent. Below 
we make the approximation that r st is constant. In reality, 
in the Epstein regime, density fluctuations will cause 
to vary with time. In the Stokes regime, both density and 
velocity fluctuations will have the same effect. However, we 
expect the essential physics to be retained if we adopt an 
appropriate time averaged value of r 8 t, so we assume this to 
have been done in what follows. 

Equation gj can be written as an equation for a 
damped harmonic oscillator randomly forced in time by the 
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part of the drag force proportional to the gaseous turbulent 
velocity field: 



dv z V z 2 u z 
— 1 h S2 z = — 

at T Bt T st 



(•>) 



The particle motion can be regarded as driven by the 
stochasting forcing term F. When the correlation time of 
the gas velocity is much smaller than the orbital period, 
this can be regarded as producing a series of weak (which 
will be the case when the stopping time is large) random im- 
pulsive changes to the velocity of the particle. FP found that 
the correlation time of the turbulent eddies is typically 0.15 
orbits, which marginally satifies this condition. Nonetheless 
we proceed with the model and later compare with simula- 
tion results. In a time interval (to,t), the change Sv z of the 
particle velocity produced by the random acceleration F can 
be written as 



Sv z = 



Fdt' = 



(6) 



Because of the stochastic nature of F, this reduces to zero 
when a suitable ensemble average is taken: {Sv z )=0. To de- 
scribe the evolution of the particle distribution, we need to 
calculate the mean square deviation of the particle velocity: 



d (Sv z y 
dt 



* u z (z,t)u z (z,t') 



dt' 



(7) 



Now, we take an ensemble average. When this is done, for 
an assumed steady state and local homogeneous turbulence, 
the result of the integral depends only on the time difference 
t — t'. Furthermore, there is no prefered time so we finally 
get 



d((Sv z ) 2 ) 
dt 



2 t* 2 f t 

= — / {u z (z,t)u z (z,0)}dt= — / S zz (t)dt. 

r st JO T st JO 

(8) 

Here S zz (t) is the velocity correlation function used in FP 
where the reader is referred for more discussion of the above 
aspects. At large times, we thus obtain the limiting form 



d((Sv z y) 2Pt(oo) 
dt ~* rl 



(9) 



where Dt(oo) = {u z (z, At)u z (z,0))d(At) defines an ef- 
fective diffusion coefficient and was already used by FP in 
a different context in their discussion of small particles well 
coupled to the gas, where it is the spatial diffusion coef- 
ficient. Integrating the last equation, we finally obtain for 
large time t 



<(*«.)"> 



2n 2L» T (oo)t 



2.3 Fokker— Planck equation 

The next step is to determine the particle distribution func- 
tion / = f(z,v z ,t), from which quantities such as the 
vertical semi-thickness can be determined. Its evolution in 
phase space can be described by the one-dimensional Fokker- 
Planck equation. For the system we study, the particles are 
subject to frictional drag and gravity due to the central star 
as written in equation while their velocities diffuse in 
velocity space with diffusion coefficient D v given by equa- 
tion IIUI I. Then, the Fokker-Planck equatio n takes the form 
<Chandrasekhadll949l I Johnson et afll2006h : 

where b = — Q, 2 z — v z /r s t is the non-stochastic acceleration. 
The derivation of equation 11 1 2> is given in Appendix 1X1 

We seek equilibrium solutions, expected to be attained 
on a time scale comparable to the stopping time, for the 
distribution function found by setting Jj = in Eq. 11211 . 
We find a solution for / of the form 



where 



(13) 



(14) 



L>t(oo) ' 

Thus the thickness of the solid body layer Hd is written as 



We note that this form of Hd as the parameter of a Gaussian 
distribution in z is of the same f rom as applies to the strongly 
coupled small particle case Csee lDubrulle et alll995l or FP). 
Thus the expression properly extends to all particle sizes. 



(10) 



3 NUMERICAL SIMULATIONS 

In this section, we present a set of numerical simulations that 
confirm the results of the discussion above. Those are local 
simulations of stratified discs that use the shearing sheet 
approximation llHawlev et al.lll995l : Istone et alJll99ch . The 
disc becomes turbulent in the presence of a weak magnetic 
field. When turbulence is well established, we follow the evo- 
lution of dust particles and larger solid bodies using two dif- 
ferent algorithms and compare the steady-state thickness of 
the solid layer to the results presented above. Below, we de- 
scribe the set-up of the simulations before going on to give 
our results. 



The linear dependence in time in this expression shows that 
the particles experience a random walk in velocity space, 
with an effective diffusion coefficient 



D v = 



Dt(oo) 



(11) 



We comment that if we set t ~ r B t in equation 1101 . we 
expect that the mean square velocity dispersion ((5v z ) 2 } ~ 
2DT(oo)/r st , a form found bv lVoelk et al.l 1^980) from con- 
siderations of stochastic forcing and a specific turbulence 
model. 



3.1 Method 



We used the ZEUS code iStone fc Normanl Il992elbh to 
evolve the standard MHD eq uations in a shearing box 
llGoldreich fc Lvnden-Beljl965l) . For the gas, we used a set- 
up identical to that described in FP. In standard Cartesian 
coordinates (x, y, z), the box has a size (H, 2nH, 6H) and we 
use a resolution (N x ,N y ,N z ) = (32, 100, 196). We comment 
that the shearing box represents a local patch of a differen- 
tially rotating disc and has been used by many authors (eg. 
lHawlev et al.lll995t IStone et al.lll996l . FP) for the purposes 



4 Carballido, Fromang & Papaloizou 




a -0 20 10 40 SO 60 /0 80 90 ! 00 ' 20 JO 40 KJ 60 /0 80 9C 100 




"<no (orcii'ls) Time (e'S>*<g) 



Figure 1. Root mean square vertical displacement as a function of time for the couplings Qr B t = 0.02, 0.2, 2, 5, 10, 30, 100 and 300 
(given at the upper left hand corner of each plot). Note the different scale used on each plot. 



of simulating MHD turbulence. The coordinate origin is at 
the centre of the box with the x and y axes pointing in the 
outward radial direction and the direction of rotation respec- 
tively. The length scale is H, the time scale is set by , 
and mass scale is determined by the gas density scale. Here 
the self-gravity of the gas is neglected. There is no informa- 
tion concerning the distance to the centre of the disc other 
than it should be much larger than H. However, the shearing 
box can be thought of as providing a detailed representation 
of a local region of a disc for which H and the distance to 
the centre are determined by global considerations. 

In our simulations, the equation of state is isothermal. 
The disc, initially in hydrostatic equilibrium in the vertical 
direction, is threaded by a magnetic field with zero net flux, 
such that the ratio of thermal to volume averaged magnetic 
pressure, (3, is equal to 400. Because of the MRI, MHD tur- 
bulence develops and is maintained for at leas t 100 orbits. 
In ag reement with previous published results llStone et alJ 
1996), we found that a, the ratio of the volume averaged to- 
tal stress (the sum of the (x,y) components of the Maxwell 
and Reynolds stress tensors) to Po, which governs the angu- 
lar nomentum transport, is typically 0.015 throughout this 
type of simulation (see FP for details). 

This model was then used as a basis to investigate the 
evolution of an ensemble of solids of different sizes. To do so, 
we used two different algorithms. First, we adopted an N- 
body approach well suited to the bodies of large size that we 
want to consider. In that case, the velocities and positions 
of individual particles are updated in two steps: the gas ve- 
locities are first interpolated at the position of the particles 
using a bilinear interpolation method. Then, a second order 



Runge-Kutta method is used to update the velocities and 
positions of the particles in time according to their equation 
of motion. For the smallest particles that we studied, we also 
used the two-fluid algorithm presented by FP, in which the 
dust component is modelled as a second pressureless fluid. 
In both cases, neither the back reaction of the particles on 
the gas nor inter-particle interactions are included. 

The N-body approach was used to investigate the evolu- 
tion of solid bodies of 8 different sizes. In each case, the par- 
ticles were assumed to be in the Epstein regime, so that r st is 
evaluated according to equation @. To label the model, we 
used the value of the dimensionless parameter f2r s t, where 
r st is evaluated in the midplane of the unperturbed disc 
(i.e. the value of the density is that given by the hydrostatic 
equilibrium) . Following this convention, the 8 models we ran 
correspond to ST2r st = 0.02, 0.2, 2, 5, 10, 30, 100 and 300. At 
t — 10 orbits, we introduce 100 particles for each of the mod- 
els described above. They are distributed in the midplane of 
the disc where they form a horizontal lattice consisting of 25 
columns and 4 rows of particles. Their initial velocities are 
set equal to zero. The system is then evolved for a further 
90 orbits. This long time integration enables a good sam- 
pling of the flow properties for each particle. This is why a 
larger number of particles is not required. Instead, we find 
that time-averaging the results over the duration of the runs 
provides a suitable statistical representation of the effect of 
MHD turbulence on the particle distribution. 

In addition, two runs were performed with the two-fluid 
method mentioned above. They are limited to f2r s t = 0.02 
and 0.2 (indeed, in order for the two-fluid approach to be 
valid, we require f2r s t < 1). As for the N-body calculation, 
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Figure 2. Vertical profile for the dust to gas ratio in the two 
models, f2r s t = 0.02 (solid line) and 0.2 (dashed line), in which 
the solid component is treated as a pressureless fluid. Both curves 
are time averaged between t = 65 and t = 85 orbits. The dotted 
line is a Gaussian fit to the data, using a thickness H d = 0.55 and 
0.18 respectively. 



Figure 3. Root-mean-square vertical displacement as a function 
of dimensionless stopping time. The filled circles (Epstein drag) 
and the asterisks (Stokes drag) represent data from simulations of 
individual particles, while the two diamonds were obtained from 
calculations regarding the dust component as a fluid. The straight 
line represents the equation (z 2 ) 1 / 2 = 0.08 (first) -1 / 2 ■ 



we assume that the Epstein law applies. The set-up we used 
for these two-fluid models is identical to that of FP: at t — 
20 orbits, the dust particles form a Gaussian thin disc such 
that Hd = 0.2H and their evolution is then followed for 80 
orbits. 



averaging the root mean square deviation of the particles be- 
tween 50 and 100 orbits. The two diamonds are the data ob- 
tained from the two-fluid calculations described above (the 
stars appearing in figure [3] are further discussed in the dis- 
cussion section below) . The straight line shows the equation 



3.2 Results 

The results of the 8 models that use the N-body approach 
are illustrated in figure Q For each of them, a curve is plot- 
ted that shows the variations of the averaged root mean 
square displacement (Z 2 } 1//2 over all the particles with time. 
The value of fir st corresponding to each case is given in 
the upper left hand corner of each plot. For convenience we 
adopt a system of units such that H = 1. However, note the 
different scales on the vertical axes of the four panels. In 
each model, (Z 2 ) 1/2 first increases before oscillating aroung 
a mean value. This mean value represents the average plan- 
etesimal disc semi-thickness. It is quickly well defined for 
the smallest particles (fir s t 30). For fir s t — 100 and 300, 
more time is required before reaching an equilibrium. This is 
because the interaction between the particles and the gas is 
very weak in that case (r st = 15 and 45 orbits respectively). 

The results of the models that use the two-fluid ap- 
proach are illustrated in figure This shows the vertical 
profile of the dust to gas mass density ratio, averaged in 
time between t = 65 and t = 85 orbits for the two models: 
fir s t = 0.02 (solid line) and fir st = 0.2 (dashed line). The 
dotted curves show the function 

with Hd/H — 0.55 and 0.18 respectively. These give the 
equilibrium dust disc semi-thickness for these two models. 

In Fig. [3] we summarize the results of all models by 
plotting the root mean square vertical extents as a function 
of fir s t. The filled circles correspond to the simulations per- 
formed with the N-body algorithm. They are obtained by 



<2 2 ) 1/2 = 0.08(fir st )- 1/2 . (17) 

which is the best fit to the data. It is seen to fit nicely 
the results of the simulations. This is expected according to 
the discussion given in section |^] Indeed, using the relation 
H — c s /fi, equation 115H can be expressed as 

From the correlation function of the turbulent velocities, FP 
found Dt(oo) = 5.5 x 10~ 3 c s H, which gives 

^ = 7.4x 10~ 2 (first)" 1/2 . (19) 
ti 

This is in very good agreement with the results of the sim- 
ulations given by equation 1171 , so validating the discussion 
presented in Section 

The evolution of (« 2 ) 1/2 , (u 2 ,) 172 and {v 2 ) 1/2 , being the 
the root mean square velocity components in the x, y and z 
directions, respectively, as a function of time is illustrated 
in figure [I] for the particles with fir s t = 10. Other values of 
the latter parameter give similar results. After initial tran- 
sients, noting that the y component is measured relative 
to the local mean shear, the ratios (v 2 ) 1 ^ 2 / (v 2 ) 1 ^ 2 ~ 0.5 
and (v 2 ) 1 / 2 / (v 2 ) 1 / 2 ~ 0.4 are established. These can be 
understood as corresponding to the root mean square tur- 
bulent gas velocity components which show similar ratios 
jHawlev et alJll995t) , In a shearing box simulation which 
reaches a quasi steady state such ratios are established and 
maintained. Restoring units, these velocity measures scale 
as Q.H. 
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Figure 4. The form of the root mean square velocity components 
divided by Q as a function of time for the particles with Qr B t = 
10. After initial transients, <f 2 ) 1/2 ~ 2.5(v 2 ) 1 / 2 and (v 2 ) 1 / 2 ~ 
1.25{u 2 ) 1 / 2 . The y component is measured relative to the local 
mean shear. 



4 DISCUSSION 

We have presented an analytical estimate of the thickness 
of the solid layer formed by sedimenting solid bodies in a 
turbulent protoplanetary disc. The stochastic dust vertical 
displacement in the presence of turbulence is modelled using 
a Fokker-Planck equation, which describes the evolution of 
the particle distribution function in phase space. The steady 
state solution of this equation provides an expression for the 
thickness Hd of the disc of solids. We found 

H d /H oc (Qr st )- 1/2 . (20) 

To confirm this analytical result, we performed simulations 
of dust particles embedded in a vertically stratified and tur- 
bulent local model of a protoplanetary disc. We found the 
numerical values of the particle height dispersions to agree 
well with the Fokker-Planck result. These simulations were 
done utilising Epstein drag. In the following, we discuss the 
expected differences that will arise when Stokes drag is used. 

4.1 The Stokes drag law 

In the Stokes regime, the stopping time r st is given by equa- 
tion For large particles, mostly confined in the vicinity 
of the equatorial plane, we should be in the regime in which 
\v\ « \u\. In that case, an averaged expression for the 
stopping time can be written as 

<T- >(»). «_ ww^y 1 , (21) 

3Gl5 pc s \ c s J 

where p is the mean value of the turbulent density field and 
(Su 2 ,) 1 ^ 2 is the root mean square turbulent gas vertical ve- 
locity fluctuation. We note that this expression gives a very 
naive estimate of the average value of (r st )' st '. Indeed, r st , 
as defined in equation depends on the turbulent veloc- 
ity. As such, the theory developed in section T2.2I should be 
modified and thus higher order velocity correlations would 
appear. Such an analysis in detail is beyond the scope of 
this paper, but the replacement of the quadratic velocity 



dependence in the drag law by the product of a linear one 
and an average in the simple estimate is likely to underes- 
timate the particle force fluctuations and accordingly un- 
derestimate the diffusion coefficient. As we will see below, 
however, such effects remain small. In equation 12 U . we have 
also considered Cd to be constant as in the large lZ e limit. 
However, a dependence of Cd on (Su 2 ) 1 ' 2 could easily be 
incorporated. 

Using equation 1211 . the dimensionless parameter 
(f2r s t) ( ' S '- ) can be written in terms of the dimensionless pa- 
rameter (f2-r s t corresponding to the Epstein regime, for 
which 

{nr st ) (Ep) = Bfl . (22) 

PCs 

Indeed, we have: 

<fiTst)<St) = 3C^ rst)(BP) (^S^-)" 1 ■ (23) 

In our numerical simulations, we found (<5u 2 ) 1/2 to be ~ 
0.13c s in th e disc midplane. This is in agreement w ith previ- 
ous results (IStone et alJll99rJ: iTurner et aT]l2006ft. Further - 
more, for the largest solid bodies. Iweidenschillind (ll9T7t> 
gives Cd = 0.44. In that case, equation I23H gives 

(fir st ) (st) ~45<fir st ) (Ep) (24) 

This relation indicates that particles in the Stokes regime 
should behave similarly to particles in the Epstein regime, 
but with a value of the dimensionless parameter first that 
is 45 times larger. 

To check this, we performed two simulations in which 
the particles follow the Stokes drag law, with Cd = 0.44, and 
we measured the equilibrium semi-thickness of the particle 
discs. The values of the dimensionless parameter (first)' Bp ' 
for these runs were taken to be equal to 0.24 and 2.4 re- 
spectively. According to the discussion above, the effective 
dimensionless parameters (fir st ) < ' St ' ) are respectively ~ 11 
and ~ 110. The measured particle disc thicknesses are rep- 
resented in figure |3] with stars. 

Although they lie close to the solid line (the vertical 
shift is between 50% and a factor of 2), there seems to be 
a systematic difference with the Epstein results. This arises 
because, as we discussed above, the expression for the dif- 
fusion coefficient should be adjusted upwards. Because of 
the quadratic dependence on velocities of the drag force, we 
expect the points derived using the Stokes law to lie system- 
atically above the equivalent Epstein points, as is observed 
in figure |3] Taking the limitations of the simple theory into 
account, the broad agreement of these simulations with the 
Epstein results is quite satisfactory and confirms our analy- 
sis. 

4.2 Gravitational instability 

Given the above results it is of interest to investigate 
whether gravitational instability of a sedimenting layer in 
a protoplanetary disc with the degree of turbulence consid- 
ered here could be a possibility. We begin by remarking that 
studies of the degree of ionisation indicate that this may be 
adequate to allow the MRI to operate both at small dis- 
tances < 0.2 AU and large distanc es > 10 — 50 AU from the 
central star jFromang et alJl200Sft . depending on the disc 
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model. In the intermediate range, a dead zone is expected 
to be present and affect the results presented below. We will 
therefore concentrate on the outer parts of the disc in the 

foll owing. 

IWeidenschillind jl994l Il995l) has discussed the possible 
role of gravitational instability in the outer solar system in 
the formation of cometary aggregates. He considered bodies 
with a ~ 10 3-4 cm with a velocity dispersion determined 
by their radial drift driven by an underlying pressure gradi- 
ent. Whether such large bodies may be for med through the 
aggreg ation of small particles is uncertain. IWeidenschillind 
( 1997) put forward a model applicable to a disc, without 
MHD turbulence, where bodies in this size range could ag- 
gregate by collisional sticking in the outer solar system on a 
10 5 year time scale. It is unclear to what extent the presence 
of MHD turbulence would extend this. The work described 
bv lWeidenschilling fc Cuzzil il993T) . appropriate to the inner 
solar system, indicates that the time scale may be extended 
by about an order of magnitude, in which case they could 
form within expected lifetimes of protostellar discs. 

Here we shall assume that bodies in the 10 4 cm size 
range can form and consider the effect of turbulence on the 
possibility of gravitational instability in the outer protoplan- 
etary disc. In the first instance we neglect radial drift as- 
suming that we are considering an ensemble of solids near a 
pressure maximum and return to this aspect later. 

For the disc model, we adopt a surface density distribu- 
tion as a function of radius r that follows that of the min- 
imum mass solar nebula 



density of the solid disc, <r p , is given by 



— !I2 
7p ~ X V5AU 



-3/2 



The surface 



(25) 



where <ro is equal to 150 g.cm -2 for the minimum mass solar 
nebula and \ is the g as 1° dust ratio. The optical depth of 
the disc of solids is 



(26) 



3ov, 

T = . 

4p s o 

For the surface density profile I25H this gives 

T= 4^l5Au) ■ (27) 

Thus for r ^ 5 AU, and a > 1 cm, the disc of solid material 
will be optically thin if we assume x > 100 (r/5AU) -3 ^ 2 . It 
will then not be able to act collectively as a fluid dust layer 
as in models that generate turbulence through interaction 
with the gas fee. iGoldreich fc Wardlll9 73: Wcideiischilline 
I1980D . 

Gravitational in stability of the particle layer can occur 
if the iToomrd lll964l) criterion is satisfied: 



c p" < 1 
nGcr D 



(28) 



where c p is the particle velocity dispersion in the radial or 
x direction and G is the gravitational constant. The ver- 
tical scale height of the particle layer is = (w 2 ) 1 ' 2 
Because the turbulent velocity fluctuations are larger in the 
radial direction than in the azimutal and vertical directions, 
the induced root mean square velocity component fluctua- 
tions induced in the solid bodies exhibit the same ordering 
(see figure Hand above). We found {vl) 1/2 ~ 2.5{v 2 ) 1/2 



and {vl) 1 ' 2 ~ 1.25<« 2 > 1/2 . Thus c p = (v 2 x )^ 2 ~ 2.5<« 2 ) 1/2 . 
Equation 1281 can then be expressed as 



Ha 3.5 x IP" 6 / r 
r < X ^5AU 



1/2 



1 g.cm 



(29) 



where the mass of the central star has been taken equal 
to one solar mass and equation 1251 has been used. Using 
Eq. I19H . which gives the height of the dust layer determined 
by turbulent diffusion, we find that 



flr st > 1.1 X 10V 



— V 

0.05r / 



— V 

5AU/ 



(To 



1 g.cm 



(30) 

Assuming the particles are in the Epstein regime (at 50 AU, 
this is valid for particles smaller than approximately 1000 
metres) together with p s = 2 g.cm -3 , a condition on their 
size can be found from equation (J5| in the form 



a > 2.8x 



H 



0.05r/ V5AU 



r \-s/a 



) 



(To 



g.cm" 



km . 



(31) 

Thus at 50 AU, for H/r = 0.05, gravitational instability 
becomes possible for a > 600 m for the minimum mass solar 
nebula if x = 100. But note that this quantity decreases 
with both (To and r and increases with \- F° r example, in 
a disc of only 3 times the minimum mass solar nebula, if 
X = 50, solid boulders larger than ~ 50 metres will become 
gravitationally unstable at 50 AU. 

The length scale associated with the instability is fc" 1 = 
Hd, with a characteristic mass M p ~ 67ro- p // 2 which gives 



Mr, 



1.3 x 10 1 



r \ »/2 
5AU7 



(To 



g.crn" 



(32) 



corresponding to a solid body of radius ~ 
5.4(r/5AU) 1/2 x _1 (cro/lg.cm- 2 ) km. At 50 AU, for 
the minimum mass solar nebula and when x — 100, this 
corresponds to a radius o f ~ 26 km, sim i lar to typical radii 
of Kuiper belt objects JLuu fc Jewitd 12002ft . In regions 
where oo/x or equivalently solid surface density is locally 
enhanced, this size is accordingly increased. For example, 
for an an enhancement of a factor of 10, the masses involved 
in the gravitational instability would correspond to bodies 
of ~ 260 km. 

As noted bv IWeidenschillind <ll995T) this form of insta- 
bility does not mean that gravitational collapse ensues but 
that a g ravitationally boun d cluster of total mass M p forms 
(see also lTanga et all2004h . This then has to evolve through 
further settling and accumulation. Nevertheless, the analysis 
presented here shows that MHD turbulence does not neces- 
sarily prevent bodies of plausibly up to few 100 km in radius 
forming by gravitational instability in the outer regions of 
protoplanetary discs. 

Note that, in drawing this conclusion, we completely 
neglected the effect of radial migration and concentrated 
only on vertical settling. This is because, in general, the ra- 
tio of the radial drift speed to orbital speed, induced by a 
radial pressure gradient, bei ng of order (H/r) 2 /(fir s t) (see 
iPapaloizou fc TerauemlEoOrl and references therein), is sig- 
nificantly less than Hd/r near marginal stability for the con- 
ditions we consider, thus it can be safely neglected. 
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APPENDIX A: 
EQUATION 



THE FOKKER-PLANCK 



The purpose of this appendix is to present a derivation of 
the Fokker-Planck equation I12H that we used in section|5|to 
estimate the pl anetesimal disc thickn ess. This derivation is 
largely based on lChandrasekhar] lll949fl in which applications 
of the theory of Brownian motion to dynamical friction and 
stellar dynamics are discussed. 

Consider the motion of a given solid body in the vertical 
direction. During a small time At , it will undergo a change 
in velocity Av z . Some of this change is due to stochastic 
forces that arise from the turbulent gas, and which have 
correlation times that are assumed to be smaller than the 
time required for typical velocities of the solid bodies to 
change significantly. Also velocity changes induced during 
a characteristic correlation time are assumed to be small. 
This is well satisfied for the large particles we consider in 
this paper. From equation @, we know that Av z relates to 
the stochastic forcing term F through 



F(At) = Av z + Q 2 zAt + (v z /r)At 



(Al) 



The probability distribution function of F(At) is known 
from the theory of Brownian motion: 

F(At) 2 ' 



^{z,v z ,Av z ,At) 



v / 47rqAt 6XP V 4(?A£ 

where F(At) is given by equation iAll and q = D v is the 
diffusion coefficient associated with the diffusive evolution of 
the velocity. We demonstrated in section T2.2l that it relates 
to the velocity correlation function through equation 1111 . 

Consider now the ensemble constituted by all the par- 
ticles in the volume of interest. Let f(z,v z ,t) be its distri- 
bution function. Under the assumption that the turbulence 
is a stationary st ochastic process wit h finite amplitude and 
correlation time djohnson et aljl2006l) . its time evolution re- 
lates to the probability function 9 through: 

/ + oo 
f(z - v z At,v z ~ Av z ,t) x 

$(z-v z At,v z - Av z ,Av z , At)d(Av z ) (A3) 
This equation can be expanded to first order in At to give 



(A2) 



Of 



dt z dz dv z 



(Av z ) al 
At 



1 d 2 



{Av 2 ) av 
At 



(A4) 



where the quantities (Av z ) av and (Av 2 ) av are defined by 
the relations: 



{Av z )av 



Av z ^(z, v z ,Av z ,At)d{Av z 



(A5) 



/ + oo 
(Av z ) 2 V(z,v z ,Av z ,At)d(Av z ). (A6) 
- oo 

Using equation 1A2L both quantities can be evaluated to 
first order in At: 



{Av z )av 

(Av 2 z )av = 2qAt 



Tst 



Q. 2 z)At, 



(A7) 
(A8) 



After eliminating At, we finally obtain equation 11 2ti by sub- 
stituting equations 1A7II and iASl in equation 1A4I . 



